# Bernd Beber, Michael J. Gilligan, Jenny Guardado, and Sabrina Karim
# R x64 3.4.4 code to replicate analysis for "The Promise and Peril of Peacekeeping Economies"
# Code depends on Stata output; run BGGK_ISQ_code.do first

# Figure 1
    # Durch's case 1 (2010), i.e. with a steady percentage of mission procurement's local content over time (actual, not survey-based)
    # display missions with complete series (i.e. not Sudan)
    # input data
    Burundi <- c(2.9, 2.5, 1.6, .74)
    DRC <- c(1.1, 1, 1, .94)
    CDI <- c(.1, .2, .2, .15)
    Haiti <- c(.7, .9, .8, .85)
    Kosovo <- c(3, 2.7, 2.1, 2.17)
    Liberia <- c(5.9, 4.7, 4.2, 3.79)
    Sierra_Leone <- c(.9, .42, .3, .24)
    Timor_Leste <- c(2, 1.08, 2.9, 3.13)
    Years <- c(2004.5, 2005.5, 2006.5, 2007.5)
    # generate plot
    pdf("Fig_1_GDP_impact.pdf", 10, 8)
        par(mar=c(5, 4, 2, 2) + 0.1)
        plot(NA, ylim=c(0,6), xlim=c(2004,2008), ylab="Local impact as percent of GDP (Durch 2010, case 1)", xlab="", bty="n")
        lines(Years, Burundi)
        lines(Years, DRC)
        lines(Years, CDI)
        lines(Years, Haiti)
        lines(Years, Kosovo)
        lines(Years, Liberia)
        lines(Years, Sierra_Leone)
        lines(Years, Timor_Leste)
        text(2004.45, Burundi[1], "Burundi", adj=c(1,.5))
        text(2004.45, DRC[1], "DR Congo", adj=c(1,.5))
        text(2004.45, CDI[1], "Cote d'Ivoire", adj=c(1,.5))
        text(2004.45, Haiti[1], "Haiti", adj=c(1,.5))
        text(2004.45, Kosovo[1], "Kosovo", adj=c(1,0))
        text(2004.45, Liberia[1], "Liberia", adj=c(1,.5))
        text(2004.45, Sierra_Leone[1], "Sierra Leone", adj=c(1,.5))
        text(2004.45, Timor_Leste[1], "Timor-Leste", adj=c(1,.5))
        lines(Years, rep(.5, length(Years)), lty=2, lwd=2)
        text(2007.55, .5, "Marshall Plan\nfor Germany\n(excl. debt\nrelief)", adj=c(0,.5))
        lines(Years, rep(5.5, length(Years)), lty=2, lwd=2)
        text(2007.55, 5.5, "2009 U.S.\nstimulus\npackage", adj=c(0,.5))
    dev.off()

# Figure 2
    # read in cross-country data
    require("readstata13")
    data_macro <- read.dta13("BGGK_ISQ_data_macro.dta")
    attach(data_macro)
    # generate plot
    pdf("Fig_2_El_Salvador.pdf", 12, 8)
        plot(1975:2005, growthpc[countryname=="El Salvador" & year>1974 & year<2006], type="n", xlab="Year", ylab="GDP per capita growth (annual %)")
        conflictyears <- (1975:2005)[as.logical(conflict[countryname=="El Salvador" & year>1974 & year<2006])]
        for(i in conflictyears) polygon(c(i-.5, i+.5, i+.5, i-.5), c(rep(par("usr")[3], 2), rep(par("usr")[4], 2)), border=NA, col="grey65")
        unyears <- (1975:2005)[as.logical(unpkm[countryname=="El Salvador" & year>1974 & year<2006])]
        for(i in unyears) polygon(c(i-.5, i+.5, i+.5, i-.5), c(rep(par("usr")[3], 2), rep(par("usr")[4], 2)), border=NA, col="grey80")
        box()
        lines(cbind(1975:2005, 0), lty=2)
        lines(1975:2005, growthpc[countryname=="El Salvador" & year>1974 & year<2006], lwd=1.5)
        text(1985, par("usr")[3] + 1, "Salvadoran Civil War")
        text(1993.5, par("usr")[3] + 1, "ONUSAL")
    dev.off()

# Figure 4
    # read in data from Liberia
    data_micro <- read.dta13("BGGK_ISQ_data_micro.dta")
    attach(data_micro)
    # generate plot
    pdf("Fig_4_share_foreign.pdf", 8, 6)
        par(mar=c(5.1, 5.1, 2, 1))
        hist(share_foreign, col="grey90", main="", xlab="Share of customers that are foreign")
    dev.off()

# Figure 5
    # occupations by composition of customer base
    m <- matrix(NA, 9, 3)
    for(i in 1:9) {
        m[i, 1] <- sum(occupation[share_foreign==1]==i, na.rm=T)
        m[i, 2] <- sum(occupation[share_foreign<1 & share_foreign>=.1]==i, na.rm=T)
        m[i, 3] <- sum(occupation[share_foreign<.1]==i, na.rm=T)
    }
    m <- apply(m, 2, function(x) x/sum(x))
    # add labels
    rownames(m) <- c("Managers", "Univ.-educ. professionals", "Other professionals", "Clerical support workers", "Service and sales workers", "Skilled agricultural workers", "Craft workers", "Machine operators", "Elementary occupations")
    colnames(m) <- c("Only foreign", "Local and at least 10% foreign", "Less than 10% foreign")
    # reverse order for display
    m <- m[nrow(m):1, ]
    # generate bar graphs
    pdf("Fig_5_occupations.pdf", width=8, height=8)
        # set margins: bottom, left, top, right
        par(mar=c(4.1, 3.1, 1.1, .2))
        barplot(m, yaxt="n", main="", xlab="Customer base")
        axis(2, seq(0, 1, .2), c("0%", "20%", "40%", "60%", "80%", "100%"), las=1, cex.axis=.9, hadj=1) 
        # range in bar plot is given by [space, (width + space) * number of columns], plus 4% on either end
        # by default, space=.2 and width=1, so for three columns we have columns centered at c(.7, 1.9, 3.1)
        # add labels
        for(i in 1:nrow(m)){
            if(!i %in% c(4, 9)) text(.7, sum(m[0:(i-1), 1]) + m[i, 1] / 2, rownames(m)[i], cex=.9, col=ifelse(i==1, "white", "black"))
            if(!i %in% 4) text(1.9, sum(m[0:(i-1), 2]) + m[i, 2] / 2, rownames(m)[i], cex=.9, col=ifelse(i==1, "white", "black"))
            if(!i %in% 4) text(3.1, sum(m[0:(i-1), 3]) + m[i, 3] / 2, rownames(m)[i], cex=.9, col=ifelse(i==1, "white", "black"))
        }
    dev.off()

# Figure 6
    # main source of funding
    f <- table(biz_funding)
    f <- f[!names(f) %in% c("Don't know", "Refuse", "")]
    f <- f[order(f, decreasing=T)]
    # collapse uncommon responses into "other"
    f <- c(f, sum(f[f<=5], f["Other, specify"]))
    names(f)[length(f)] <- "Other"
    f <- f[!f<=5]
    f <- f[!names(f)=="Other, specify"]
    # calculate proportions
    f <- f/sum(f)
    # adjust labels
    names(f)[names(f)=="Gift from family"] <- "Family gift"
    names(f)[names(f)=="Loan from bank"] <- "Bank loan"
    names(f)[names(f)=="Loan from savings club"] <- "Savings club loan"
    # create matrix and add column label
    f <- matrix(f, dimnames=list(names(f), "Main funding source to start business"))

    # main problem in running business
    p <- table(biz_problem)
    p <- p[!names(p) %in% c("Don't know", "No major problem", "")]
    # collapse crime and corruption
    p <- c(p, sum(p["Corruption of government officials"], p["Crime (theft, rackets)"]))
    names(p)[length(p)] <- "Crime/corruption"
    p <- p[!names(p)=="Corruption of government officials"]
    p <- p[!names(p)=="Crime (theft, rackets)"]
    # order
    p <- p[order(p, decreasing=T)]
    # collapse uncommon responses into "other"
    p <- c(p, sum(p[p<=5], p["Other, specify"]))
    names(p)[length(p)] <- "Other"
    p <- p[!p<=5]
    p <- p[!names(p)=="Other, specify"]
    # calculate proportions
    p <- p/sum(p)
    # adjust labels
    names(p)[names(p)=="Problems with equipment or spare parts"] <- "Problems with equipment"
    # create matrix and add column label
    p <- matrix(p, dimnames=list(names(p), "Most serious problem in running business"))

    # generate bar graphs
    pdf("Fig_6_business.pdf", width=10, height=8)
        par(mfrow=c(1, 2))
        # set margins: bottom, left, top, right
        par(mar=c(3, 3.1, 2, .2))
        barplot(f, yaxt="n")
        axis(2, seq(0, 1, .2), c("0%", "20%", "40%", "60%", "80%", "100%"), las=1, cex.axis=.9, hadj=1)
        for(i in 1:nrow(f)){
            text(.7, sum(f[0:(i-1), 1]) + f[i, 1] / 2, rownames(f)[i], cex=.9, col=ifelse(i %in% c(1, 2), "white", "black"))
        }
        par(mar=c(3, 2, 2, .2))
        barplot(p, yaxt="n")
        for(i in 1:nrow(p)){
            text(.7, sum(p[0:(i-1), 1]) + p[i, 1] / 2, rownames(p)[i], cex=.9, col=ifelse(i %in% c(1, 2), "white", "black"))
        }
    dev.off()

# Figure 10
    # UNMIL popularity
    # read estimation output
    out <- read.csv("_out.csv")
    attach(out)
    # generate graph
    pdf("Fig_10_popularity.pdf", 12, 7)
        par(mgp=c(3,2,0))
        # for regression estimates, see do-file
        x <- barplot(p1, names.arg=c("Wants UNMIL to stay", "UNMIL has improved\ncity's security", "UNMIL has\nimproved own security", "Positive non-security\nconsequences", "No negative\nconsequences"), col=NA, density=20, ylim=c(0, 1))
        barplot(p0, add=T, col="white")
        text(x, .04 + p1, paste("+", round(100 * (p1 - p0), 1), "%"))
    dev.off()
